using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;

namespace SmartMathLibrary.QuadPack
{
    /// <summary>
    /// Integration over finite intervals. Accepts a list of known singularities.
    /// </summary>
    [Serializable]
    public static class DqagpClass
    {
        /// <summary>
        /// Adaptive integration routine which handles functions to be integrated between two finite bounds. 
        /// The adaptive strategy compares results of integration over the given interval with the sum of results 
        /// obtained from integration over a bisected interval. Since error estimates are available from each 
        /// regional integration, the region with the largest error is bisected and new results are computed. 
        /// This bisection process is continued until the error is less than the prescribed limit or convergence 
        /// failure is determined.
        /// </summary>
        /// <param name="f">Double precision function to be integrated.</param>
        /// <param name="a">Lower limit of integration.</param>
        /// <param name="b">Upper limit of integration.</param>
        /// <param name="npts2">The number equal to 2 more than the number of sinularities.</param>
        /// <param name="points">The vector of dimension npts2, the first (npts2-2) elements of which are the user provided interior break points.</param>
        /// <param name="epsabs">The absolute accuracy requested.</param>
        /// <param name="epsrel">The relative accuracy requested.</param>
        /// <param name="abserr">The requested absolute error of the result.</param>
        /// <param name="neval">The needed function evalations.</param>
        /// <param name="ier">The requested errorcode.</param>
        /// <returns>The result of the numerical integration.</returns>
        public static double Dqagp(HardRealFunction f, double a, double b, int npts2, double[] points, double epsabs,
                                   double epsrel, ref double abserr, ref int neval, ref int ier)
        {
            double abseps = 0;
            double[] alist = new double[QuadPackConstants.LIMIT];
            double area;
            double area1;
            double area12;
            double area2;
            double a1;
            double a2;
            double[] blist = new double[QuadPackConstants.LIMIT];
            double b1;
            double b2;
            double correc = 0;
            double defabs = 0;
            double defab1 = 0;
            double defab2 = 0;
            double dres;
            double[] elist = new double[QuadPackConstants.LIMIT];
            double erlarg;
            double erlast;
            double errbnd;
            double errmax;
            double error1 = 0;
            double error2 = 0;
            double erro12;
            double errsum;
            double ertest;
            double[] ndin = new double[40];
            double[] pts = new double[40];
            double resa = 0;
            double resabs;
            double reseps;
            double result;
            double[] res3la = new double[3];
            double[] rlist = new double[QuadPackConstants.LIMIT];
            double[] rlist2 = new double[52];
            double sign;
            double temp;

            int i;
            int id;
            int ierro;
            int ind1;
            int ind2;
            int ip1;
            int[] iord = new int[QuadPackConstants.LIMIT];
            int iroff1;
            int iroff2;
            int iroff3;
            int j;
            int jlow;
            int jupbnd;
            int k = 0;
            int ksgn;
            int ktmin;
            int last;
            int levcur;
            int[] level = new int[QuadPackConstants.LIMIT];
            int levmax;
            int maxerr;
            int nint;
            int nintp1;
            int npts;
            int nres;
            int nrmax;
            int numrl2;
            int limit;
            int extrap;
            int noext;

            limit = QuadPackConstants.LIMIT - 1;

            // Test validity of parameters. 	
            ier = 0;
            neval = 0;
            last = 0;
            result = 0.0;
            abserr = 0.0;
            alist[0] = a;
            blist[0] = b;
            rlist[0] = 0.0;
            elist[0] = 0.0;
            iord[0] = 0;
            level[0] = 0;
            npts = npts2 - 2;
            if ((npts2 < 2) || (limit < npts) || ((epsabs < 0.0) && (epsrel < 0.0)))
            {
                ier = 6;
            }
            if (ier == 6)
            {
                goto _999;
            }

            // If any break points are provided, sort them into an ascending sequence. 
            sign = (a < b) ? 1.0 : -1.0;
            pts[0] = (((a) < (b)) ? (a) : (b));
            if (npts == 0)
            {
                goto _15;
            }
            for (i = 0; i < npts; i++)
            {
                pts[i + 1] = points[i];
            }
            _15:
            pts[npts + 1] = (((a) > (b)) ? (a) : (b));
            nint = npts + 1;
            a1 = pts[0];
            if (npts == 0)
            {
                goto _40;
            }
            nintp1 = nint + 1;
            for (i = 0; i < nint; i++)
            {
                ip1 = i + 1;
                for (j = ip1; j < nintp1; j++)
                {
                    if (pts[i] <= pts[j])
                    {
                        goto _20;
                    }
                    temp = pts[i];
                    pts[i] = pts[j];
                    pts[j] = temp;
                    _20:
                    ;
                }
            }
            if ((pts[0] != (((a) < (b)) ? (a) : (b))) || (pts[nintp1 - 1] != (((a) > (b)) ? (a) : (b))))
            {
                ier = 6;
            }
            if (ier == 6)
            {
                goto _999;
            }

            // Compute first integral and error approximations. 
            _40:
            resabs = 0.0;
            for (i = 0; i < nint; i++)
            {
                b1 = pts[i + 1];
                area1 = Dqk21Class.G_K21(f, a1, b1, ref error1, ref defabs, ref resa);
                abserr = abserr + error1;
                result = result + area1;
                ndin[i] = 0;
                if ((error1 == resa) && (error1 != 0.0))
                {
                    ndin[i] = 1;
                }
                resabs += defabs;
                level[i] = 0;
                elist[i] = error1;
                alist[i] = a1;
                blist[i] = b1;
                rlist[i] = area1;
                iord[i] = i;
                a1 = b1;
            }
            errsum = 0.0;
            for (i = 0; i < nint; i++)
            {
                if (ndin[i] == 1)
                {
                    elist[i] = abserr;
                }
                errsum += elist[i];
            }

            // Test on accuracy. 
            //      last = nint; 
            neval = 21 * nint;
            dres = Math.Abs(result);
            errbnd = (((epsabs) > (epsrel * dres)) ? (epsabs) : (epsrel * dres));
            if ((abserr <= 100.0 * Double.Epsilon * resabs) && (abserr > errbnd))
            {
                ier = 2;
            }
            if (nint == 0)
            {
                goto _80;
            }
            for (i = 0; i < npts; i++)
            {
                jlow = i + 1;
                ind1 = iord[i];
                for (j = jlow; j < nint; j++)
                {
                    ind2 = iord[j];
                    if (elist[ind1] > elist[ind2])
                    {
                        goto _60; // use continue after debugging 
                    }
                    ind1 = ind2;
                    k = j;
                    _60:
                    ;
                }
                if (ind1 == iord[i])
                {
                    goto _70;
                }
                iord[k] = iord[i];
                iord[i] = ind1;
                _70:
                ;
            }
            if (limit < npts2)
            {
                ier = 1;
            }
            _80:
            if ((ier != 0) || (abserr <= errbnd))
            {
                goto _999;
            }

            // Initialization. 
            res3la[0] = 0.0;
            res3la[1] = 0.0;
            res3la[2] = 0.0;
            rlist2[0] = result;
            maxerr = iord[0];
            errmax = elist[maxerr];
            area = result;
            nrmax = 0;
            nrmax = 0;
            nres = -1; // nres = 0 
            numrl2 = 0; // numrl2 = 1 
            ktmin = 0;
            extrap = QuadPackConstants.FALSE;
            noext = QuadPackConstants.FALSE;
            erlarg = errsum;
            ertest = errbnd;
            levmax = 1;
            iroff1 = 0;
            iroff2 = 0;
            iroff3 = 0;
            ierro = 0;
            abserr = Double.MaxValue;
            ksgn = -1;
            if (dres > (1.0 - 50.0 * Double.Epsilon) * resabs)
            {
                ksgn = 1;
            }

            // Main loop. 
            for (last = npts2; last <= limit; last++)
            {
                // Bisect the interval with the nrmax-th largest error estimate. 
                levcur = level[maxerr] + 1;
                a1 = alist[maxerr];
                b1 = 0.5 * (alist[maxerr] + blist[maxerr]);
                a2 = b1;
                b2 = blist[maxerr];
                erlast = errmax;
                area1 = Dqk21Class.G_K21(f, a1, b1, ref error1, ref resa, ref defab1);
                area2 = Dqk21Class.G_K21(f, a2, b2, ref error2, ref resa, ref defab2);
                // Improve previous approximations to integral and error
                //  	and test for accuracy. 
                neval += 42;
                area12 = area1 + area2;
                erro12 = error1 + error2;
                errsum = errsum + erro12 - errmax;
                area = area + area12 - rlist[maxerr];
                if ((defab1 == error1) || (defab2 == error2))
                {
                    goto _95;
                }
                if ((Math.Abs(rlist[maxerr] - area12) > 1.0e-5 * Math.Abs(area12)) || (erro12 < .99 * errmax))
                {
                    goto _90;
                }
                if (extrap != 0)
                {
                    iroff2++;
                }
                else
                {
                    iroff1++;
                }
                _90:
                if ((last > 9) && (erro12 > errmax)) // last > 10 
                {
                    iroff3++;
                }
                _95:
                level[maxerr] = levcur;
                level[last] = levcur;
                rlist[maxerr] = area1;
                rlist[last] = area2;
                errbnd = (((epsabs) > (epsrel * Math.Abs(area))) ? (epsabs) : (epsrel * Math.Abs(area)));

                // Test for roundoff error and eventually set error flag. 
                if (((iroff1 + iroff2) >= 10) || (iroff3 >= 20))
                {
                    ier = 2;
                }
                if (iroff2 > 5)
                {
                    ier = 3;
                }

                // Set error flag in the case that the number of subintervals
                //	equals limit. 
                if (last == limit) // last == limit 
                {
                    ier = 1;
                }

                // Set error flag in the case of bad integrand behavior at some
                //	points in the integration range. 
                if ((((Math.Abs(a1)) > (Math.Abs(b2))) ? (Math.Abs(a1)) : (Math.Abs(b2))) <=
                    (1.0 + 1000.0 * Double.Epsilon) * (Math.Abs(a2) + 1000.0 * Double.MinValue))
                {
                    ier = 4;
                }

                // Append the newly-created intervals to the list. 
                if (error2 > error1)
                {
                    goto _100;
                }
                alist[last] = a2;
                blist[maxerr] = b1;
                blist[last] = b2;
                elist[maxerr] = error1;
                elist[last] = error2;
                goto _110;
                _100:
                alist[maxerr] = a2;
                alist[last] = a1;
                blist[last] = b1;
                rlist[maxerr] = area2;
                rlist[last] = area1;
                elist[maxerr] = error2;
                elist[last] = error1;

                // Call dqsort to maintain the descending ordering in the list of error
                //	estimates and select the subinterval with nrmax-th largest
                //	error estimate (to be bisected next). 
                _110:
                DqsortClass.Dqsort(limit, last, ref maxerr, ref errmax, elist, iord, ref nrmax);
                if (errsum <= errbnd)
                {
                    goto _190;
                }
                if (ier != 0)
                {
                    goto _170;
                }
                if (noext != 0)
                {
                    goto _160;
                }
                erlarg -= erlast;
                if (levcur + 1 <= levmax)
                {
                    erlarg += erro12;
                }
                if (extrap != 0)
                {
                    goto _120;
                }

                // Test whether the interval to be bisected next is the smallest interval. 
                if ((level[maxerr] + 1) <= levmax)
                {
                    goto _160;
                }
                extrap = QuadPackConstants.TRUE;
                nrmax = 1; // nrmax = 2 
                _120:
                if ((ierro == 3) || (erlarg <= ertest))
                {
                    goto _140;
                }

                // The smallest interval has the largest error. Before bisecting, decrease
                //    the sum of the errors over the larger intervals (erlarg) and
                //		perform extrapolation.) 
                id = nrmax;
                jupbnd = last;
                if (last > (2 + limit / 2))
                {
                    jupbnd = limit + 3 - last;
                }
                for (k = id; k <= jupbnd; k++)
                {
                    maxerr = iord[nrmax];
                    errmax = elist[maxerr];
                    if (level[maxerr] + 1 <= levmax)
                    {
                        goto _160;
                    }
                    nrmax++;
                }

                // Perform extrapolation. 
                _140:
                numrl2++;
                rlist2[numrl2] = area;
                if (numrl2 <= 1)
                {
                    goto _155;
                }
                reseps = DqextClass.Dqext(ref numrl2, rlist2, ref abseps, res3la, ref nres);
                ktmin++;
                if ((ktmin > 5) && (abserr < 1.0e-3 * errsum))
                {
                    ier = 5;
                }
                if (abseps >= abserr)
                {
                    goto _150;
                }
                ktmin = 0;
                abserr = abseps;
                result = reseps;
                correc = erlarg;
                ertest = (((epsabs) > (epsrel * Math.Abs(reseps))) ? (epsabs) : (epsrel * Math.Abs(reseps)));
                if (abserr <= ertest)
                {
                    goto _170;
                }

                // Prepare bisection of the smallest interval. 
                _150:
                if (numrl2 == 0)
                {
                    noext = QuadPackConstants.TRUE;
                }
                if (ier == 5)
                {
                    goto _170;
                }
                _155:
                maxerr = iord[0];
                errmax = elist[maxerr];
                nrmax = 0;
                extrap = QuadPackConstants.FALSE;
                levmax += 1;
                erlarg = errsum;
                _160:
                ;
            }
            ;
            _170:
            if (abserr == Double.MaxValue)
            {
                goto _190;
            }
            if ((ier + ierro) == 0)
            {
                goto _180;
            }
            if (ierro == 3)
            {
                abserr += correc;
            }
            if (ier == 0)
            {
                ier = 3;
            }
            if ((result != 0.0) && (area != 0.0))
            {
                goto _175;
            }
            if (abserr > errsum)
            {
                goto _190;
            }
            if (area == 0.0)
            {
                goto _210;
            }
            goto _180;
            _175:
            if (abserr / Math.Abs(result) > errsum / Math.Abs(area))
            {
                goto _190;
            }

            // Test on divergence. 
            _180:
            if ((ksgn == -1) &&
                ((((Math.Abs(result)) > (Math.Abs(area))) ? (Math.Abs(result)) : (Math.Abs(area))) <= defabs * .01))
            {
                goto _210;
            }
            if ((0.01 > result / area) || (result / area > 100.0) || (errsum > Math.Abs(area)))
            {
                ier = 6;
            }
            goto _210;

            // Compute global integral. 
            _190:
            result = 0.0;
            for (k = 0; k <= last; k++)
            {
                result += rlist[k];
            }
            abserr = errsum;
            _210:
            if (ier > 2)
            {
                (ier)--;
            }
            result = result * sign;
            _999:
            return result;
        }
    }
}